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Abstract — We consider the stochastic patterns of a system of 
communicating, or coupled, self-propelled particles in the pres- 
ence of noise and communication time delay. For sufficiently 
large environmental noise, there exists a transition between a 
translating state and a rotating state with stationary center 
of mass. Time delayed communication creates a bifurcation 
pattern dependent on the coupling amplitude between particles. 
Using a mean field model in the large number limit, we 
show how the complete bifurcation unfolds in the presence of 
communication delay and coupling amplitude. Relative to the 
center of mass, the patterns can then be described as transitions 
between translation, rotation about a stationary point, or a 
rotating swarm, where the center of mass undergoes a Hopf 
bifurcation from steady state to a limit cycle. Examples of some 
of the stochastic patterns will be given for large numbers of 
particles. 

I. INTRODUCTION 

The collective motion of interacting multi-particle systems 
has been the subject of many recent experimental and model- 
ing studies. It is especially astounding that numerous coher- 
ent states of great complexity can arise spontaneously in spite 
of the absence of a particle acting as a leader. The study of 
these swarming systems has proven useful in understanding 
the spatio-temporal patterns formed by bacterial colonies, 
fish, birds, locusts, ants, pedestrians, etc. [1], [2], [3], [4], 
[5], [6]. Moreover, these studies have provided valuable 
information that may be exploited in the design of systems 
of autonomous, inter-communicating robotic systems [7], [8], 
[9]. 

Investigators have used various mathematical approaches 
to study swarm systems. Some studies have preserved the in- 
dividual character of each agent in the system, using ordinary 
or delay differential equations (ODEs/DDEs) to describe 
their trajectories [10], [11], [12], [8]. Other researchers have 
proposed continuum models written in terms of averaged 
velocity and particle density fields that satisfy partial differ- 
ential equations (PDEs) [2], [3], [5], [6]. In addition, authors 
have also studied the effects of noise in the swarms and 
shown the existence of noise-induced transitions [13], [14]. 

More recently, authors have begun to study the effects 
of communication time-delays between particles. Time-delay 
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models are common in many areas of mathematical biology 
including population dynamics, neural networks, blood cell 
maturation, virus dynamics and genetic networks [15], [16], 
[17], [18], [19], [20], [21], [22], [23]. In the context of 
swarming particles, it has been shown that the introduction of 
a communication time-delay may induce transitions between 
different coherent states [14]. The type of transition is 
dependent on the coupling strength between particles and 
the noise intensity. 

Here we make a more detailed study of the bifurcation 
structure of the mean field approximation to the delay- 
coupled model proposed studied in [14] and investigate how 
the bifurcations in the system are modified in the presence 
of noise. 

II. The Swarm Model 

We consider a two-dimensional swarm with N self- 
propelling particles that are mutually attracted in a symmetric 
fashion. Additionally, we consider the case in which parti- 
cles communicate with each other with a time delay. The 
swarm is governed by the following system of ODEs: 
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for i = 1,2...,N, Here and Vj represent the position 
and velocity of the i-th particle, respectively; the strength of 
the attraction is measured by the coupling constant a and the 
time delay is uniform and given by r. The self- propulsion 
and frictional drag on each particle is given by the term 
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In the absence of coupling, particles tend to 



move on a straight line with unit speed |vj| = 1 as time goes 
to infinity. The term rj^t) = (j^- , rf^) is a two-dimensional 
vector of stochastic white noise with intensity equal to D and 
correlation functions (r]f \t)) = and (rjf\t)ri^\t')) = 
2DS(t - t / )S ll S lk for i, j = 1,2, . . .JV and i,k= 1,2. 

The coupling between particles arises from a time-delayed, 
spring-like potential. Hence, our equations of motion may 
be considered to be the first term in a Taylor expansion of 
other more general time-delayed potential functions about 
an equilibrium point. The model described by Eqs. (TTal - 
([Tbl with t = (i.e. no time delay) possesses a noise- 
induced transition whereby a large enough noise intensity 
causes a translating swarm of individuals to transition to 
a rotating swarm with a stationary center of mass [24], 
[14]. Regardless of which state the swarm is in (translating 
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Fig. 1 . Snapshots of a swarm taken at (a) t = 50, (b) t = 60, (c) t = 62, 
(d) t = 64, (e) t = 66, (f) t = 68, (g) t = 70, (h) i = 72, (i) t = 74, 
and (j) * = 76, with a = 4, N = 300, and D = 0.08. The swarm was 
in a rotational state when the time delay of r = 1 was switched on at 
t = 40. For a movie, see the relevant mpeg video. Figure reproduced with 
permission from [14]. 



or rotating), the addition of a communication time delay 
leads to another type of transition. This transition occurs if 
the coupling parameter a, is large enough. As an example, 
we consider a swarm that has already undergone a noise- 
induced transition to a rotational state before switching on 
the communication time delay. 

Figures [Ha)-[T]j) show snapshots of a swarm at t = 50, 
t = 60, t = 62, t = 64, t = 66, t = 68, t = 70, t = 72, 
t = 74, and t = 76 respectively. For these simulations, N — 
300, t = 1, D = 0.08, the noise was switched on at t = 
10 (causing the swarm to transition to a stationary, rotating 
state), and once in this rotating state, the time delay was 
switched on at t = 40. One can see that with the evolution 
of time, the individual particles become aligned with one 



another and the swarm becomes more compact. Additionally, 
the swarm is no longer stationary, but has begun to oscillate 
[Figs.Eg)-lIlj)]. 

This compact, oscillating aligned swarm state looks similar 
to a single "clump" that is described in [25]. However, where 
each "clump" of [25] contains only some of the total number 
of swarming particles, our swarm contains every particle. 
Additionally, while a deterministic model along with global 
coupling is used to attain the "clumps" of [25], our oscillating 
aligned swarm is attained with the use of noise and a time 
delay. 

III. Mean Field Approximation 

As we have shown, once the stochastic swarm is in 
the stationary, rotating state, the addition of a time delay 
induces an instability. We investigate the stability of the 
swarm by deriving the mean field equations and performing 
a bifurcation analysis. 

We carry out a mean field approximation of the swarming 
system by switching to particle coordinates relative to the 
center of mass and disregarding the noise terms. The center 
of mass of the swarming system is given by 
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We decompose the position of each particle into 
Ti = R + Sri, 



i = l,2...,N, 



where 
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Inserting Eq. (O into the second order system equivalent to 
Eqs. ([Tal)-(fTb1) with D = and simplifying one obtains 

R + Sh = (l - |R| 2 - 2R • 6±i - l^l 2 ) (R + Sr z ) 
a(N - 1 
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where we used Eq. (HJl in the form 
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Summing Eq. (O over i and using Eq. dU, one arrives at 
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Inserting Eq. (Q into Eq. (O the following equation for 
Sri is obtained: 
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for i = 1,2..., A". 

Equations (0 and ([8]) are fully equivalent to Eqs. dTab - 
(flbl when D = 0, and simply consist of rewriting the 
original system using the relationship between the particle 
coordinates r,, the center of mass R, and the coordinates 
relative to the center of mass r5r.j. This mapping has trans- 
formed the original 2N differential equations into 2N + 2 
differential equations. There is, however, no inconsistency 
since in the transformed set of equations only 2N of them 
are independent, because of the relation seen in Eq. (0J. 

We then obtain a mean field approximation by neglecting 
the fluctuation of the swarm particles, <5rj's, from the center 
of mass: 



R = (1 - |R| 2 ) R-a(R(t) -R(t- r)) , 
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where we made the approximation a 
consider the thermodynamic limit. 



N-l 
N 



a since we 



IV. Bifurcations in the Mean Field Equation 

The behavior of the system in the mean field approxima- 
tion in different regions of parameter space may be better 
understood by using bifurcation analysis. Equation ((9) may 
be written in component form using R = (X, Y) and 

R = (U, V) as 

X = U, (10a) 

U= (l-U 2 -V 2 )U - a(X - X(t-T)), (10b) 

Y = V, (10c) 

V = (1 - U 2 - V 2 )V - a(Y-Y(t -r)). (lOd) 

For all values of a and r, Eqs. dlOab - dlOdb have translation- 
ally invariant stationary solutions 



X = X , U = 0, Y = Y , V = 0, 



(11) 



with two free parameters Xq and Yq. They also have a three 
parameter family of uniformly translating solutions 



X = U t + X 



which requires 



U = U (h Y = V t + Y , V = V , 

(12) 
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and thus exists only for ar < 1. In the two-parameter space 
(a, r), the hyperbola ar = 1 is in fact a pitchfork bifurcation 



line on which the uniformly translating states are born from 
the stationary state (Xq,0,Yq,0). The other branch of the 
pitchfork is an unphysical solution with negative speed. 

Linearizing Eqs. (llOab - dTOdt about the stationary state, we 
obtain the characteristic equation 



(a(l — e~ Ar ) — A + A 2 ) = 0. 
It suffices to study the zeros of the function 



V{\) = a(l - e~ XT ) - A + X z = 0, 
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since the eigenvalues of Eqs. ( I lOab - dTOdb are obtained by 
duplicating those of Eq. (15[ . 

We now search for Hopf bifurcations in the two parameter 
space (a, r) by letting A = iui in Eq. ( fl5l l. Substitution leads 
to 



a — uj 2 — iuj = ae ZUJT . 



(16) 



Taking the modulus of Eq. ( TToT l. we find a at the Hopf point, 
an, is given by 

a 2 H = (a H -LU 2 ) 2 +io 2 , (17) 
or, considering w^0, 
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We eliminate a in Eq. ( TToT l by using Eq. ([T81 and taking 
the complex conjugate to obtain an equation for r at the 
Hopf point 

1 - uj 2 2lu 

T—2+^T—^ = ^. (19) 

1 + U) Z 1 + LO £ 

We obtain r by equating the arguments of both sides, being 
careful to use the branch of tan (9 in (0, 7r) since the left 
hand side of Eq. dl9l > is on the upper complex plane for 
uj > 0. The result is a family of Hopf bifurcation curves 
parameterized by uj: 

1+uj 2 
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These curves are shown in Fig. |2(a)| We may eliminate the 
parameter uj between these two equations and obtain 



THn(a) = 
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In spite of their appearance, the Hopf curves in Eqs. (120at - 
(I20bb and (l2Tb are in fact continuous at u> = 1 and a = 1, 
respectively (with the correct branch of tan 9 in (0, 7r)). From 
Eq. (I20at - (l20bt . we see that the Hopf frequency depends only 
on the value of a for all members in the family; it has the 
value one at a = 1 and the frequency tends to infinity as a 
grows. Interestingly, only the first Hopf curve of the family in 





Fig. 2. (a) Hopf (blue) and pitchfork (red) bifurcation curves in a and r 
space, (b) A zoom-in of the branches in the first panel displaying also the 
saddle to node transition (dashed black); the number in each region indicates 
the number of eigenvalues with a real part greater than zero with the solid 
lines as boundaries. (Color online.) 



Eq. (f2TT > is defined at a = 1/2; it has the value tho\ 0=1/2 = 
2. The point (a = 1/2, r = 2) which lies both on the first 
member of the family of Hopf curves and on the pitchfork 
branch is in fact a Bogdanov-Takens (BT) point [26], where 
uj = 0. None of the other Hopf branches meet the pitchfork 
bifurcation line since they tend asymptotically to infinity as 
a -> 1/2. 

We used a numerical continuation method (DDE- 
Biftool) [27] to calculate the pitchfork and Hopf branches 
in the (a, r) parameter space; these results are in per- 
fect agreement with our analytical calculations (results not 
shown). These numerical studies reveal that the number of 
eigenvalues with real part greater than zero is as indicated in 
Fig. |2(b)| In addition, our numerical continuation analyses 
also reveal node to focus transitions of the steady state. 
These occur at points where there are two real and equal 
eigenvalues, i.e. where T>(\) = and T>'(\) = 0, for A real. 
From V(\) = we obtain e~ rA = 1 ~^ A , which we can 
insert into T>(\) — to obtain 
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Fig. 3. Location of the dominating eigenvalues around the Bogdanov- 
Takens point at a = 1/2, r = 2. Parameter values are (a) a = 0.60, 
t = 2.0, (b) a = 0.48, r = 2.09, (c) a = 0.40, r = 2.01, (d) a = 0.53, 
t = 1.90, and (e) a = 0.55, r = 1.91. 



this gives the curve where the node-focus transitions occur: 

1 



(23) 



roots to be repeated, we set the discriminant to zero and 



V^l/4 

Moreover, from the solutions to Eq. (l2Zt we see that the 
repeated eigenvalues have positive real parts if r > 2 and 
negative real parts if r < 2. In Fig. |2(b)| we show the 
pitchfork and Hopf bifurcation curves overlaid with the node- 
focus transition curve given by Eq. (1231 , 

As seen in Fig. |2(b)| the pitchfork and Hopf branches, 
together with the node-focus transition curves split the area 
around the BT point into five different regions. The behavior 
of the leading eigenvalues (excluding the one at the origin) as 
one probes these five regions is shown in Figs. |3(a)|3(e)1 At 
a point directly to the right of the BT point in (a, r) space, 
the stationary solution has a pair of eigenvalues with positive 
real parts and non-zero imaginary parts [Fig. |3(a)] . Moving 
counter-clockwise in the (a, r) plane, the eigenvalue pair 
collapses onto the real line after crossing the upper branch 
of the node-focus transition [Fig. |3(b)| . Still moving in the 
same direction in parameter space, we observe two different 
instances of eigenvalues crossing the origin: first, the smaller 
of the two purely real and positive eigenvalues does so on 
the upper part of the pitchfork bifurcation line [Fig. |3(c)| 
and then the remaining purely real and positive eigenvalue 
crosses the origin on the lower part of the bifurcation line 
[Fig. |3(d)) . Finally, at the node-focus transition line, the two 
purely real and negative eigenvalues coincide on the real axis 
and acquire non-zero imaginary parts [Fig. [3(e)] , Continuing 
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Fig. 4. The limit cycle of the center of mass is shown through a comparison 
of analytical (solid line) and numerical ("cross" markers) values of an and 
lo for several choices of r. The analytical result is found using Eqs. J20at - 
(20b), while the numerical result is found using a continuation method [27] 
for Eq. l|9j- The inset shows the stochastic trajectory of the center of mass 
of the swarm from t = 45 to t = 90. Figure reproduced with permission 
from [14]. 



upwards in parameter space, the complex pair crosses the 
imaginary axis on the Hopf bifurcation curve, giving birth 
to a stable limit cycle. 

A. Numerical simulations 

Figure [4] shows an excellent comparison of the analytical 
result given by Eqs. ( I20ab -( l20bb with a numerical result 
which was found using a continuation method [27] for the 
mean field model for several choices of r. Furthermore, for 
t = 1, the value of coupling a at the bifurcation point 
is an ~ 3.2. This value of an corresponds very well to 
the change in behavior of the stochastic swarm (results not 
shown). 

More evidence of the Hopf bifurcation is seen in the inset 
of Fig. |4] The inset shows the stochastic trajectory of the 
center of mass of the swarm from t — 45 to t = 90 
for the example shown in Fig. Q] Once the time delay is 
switched on at t = 40 (with the swarm located at the 
center of the inset figure), the swarm begins to oscillate. 
The swarm moves along an elliptical path [the position of 
its center of mass is denoted at several times that correspond 
to Figs, [jib), Hid), [Qf), Eh), and [Qj)], until it eventually 
converges to the circular limit cycle. 

Figures |5(a)| and |5(b)| show a time series simulation of a 
swarm with N = 75 particles. Figure |5(a)1 shows the position 
components, while Fig. |5(b)| shows the velocity components., 
One can see that the swarm follows a circular-like path over 
time. A perturbation that is applied at t = 20 shows that for 
the chosen parameters, the pattern is stable in the presence 
of noise. 

The presence of noise introduces interesting switching 
behavior that make the initial conditions of the swarm critical 
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Fig. 5. The limit cycle of the center of mass is shown through the (a) 
position and (b) velocity time series of the swarm using N = 75 particles, 
a = 0.7, t = 2.2, and noise intensity D = 0.045. A velocity perturbation 
is applied at t = 20. 
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Fig. 6. Long time behavior of the mean particle alignment (defined in the 
text) for different values of noise intensity (D = <r 2 /2) and two different 
initial conditions. In panel (a), all particles start off from the origin with 
equal velocity vectors; in panel (b), all particles start from rest, distributed 
uniformly over the unit square. For these simulations, N = 150, a = 2 
and t = 2. The time-delay is turned on at t = 50, and the simulations ran 
until t = 300. 



in determining the long time behavior of the system. To 
demonstrate this, we have performed a series of simula- 
tions for different noise intensities and two different initial 
conditions: (i) all particles start at the origin with unit x 
and y speeds [Fig. |6(a)| and (ii) all particles are distributed 
uniformly over the unit square and start from rest [Fig. |6(b)| . 
The simulations are run until t — 300 using a coupling 
constant a = 2 and a time-delay r = 2 which is turned 
on at t = 50. Our simulations reveal that in the long time 
limit and for small values of noise, the swarm converges 
to either a compact state that rotates as a whole [case (i)] 
or to a ring state with particles going both clockwise and 
counterclockwise [case (ii)]. The asymptotic behavior of 
the system is readily visualized by calculating the mean 
alignment of the swarm particles. We quantify this mean 
alignment of the swarm by calculating the cosine between 
the directions of the i-th particle and the center of mass, 
cos$i = • R)/(|rj||R|), and then averaging over all 
particles and over the last 100 time units of simulation. 
Figure |6(a)| shows that in case (i) the particles converge 
to the compact, aligned state for low and moderate noise 
intensities. However, this state is broken up at high noise 
levels (cr w 0.8). In contrast, Fig. |6(b)| shows that in case 
(ii) the particles converge to a ring for small values of noise 
(a < 0.25), evidenced by the low values of the mean particle 
alignment in Fig. |6(b)| but converge to the aligned case 
for higher values of noise (a > 0.25). Observing the full 
simulation runs in detail (not shown) reveals a switching 



behavior: for case (ii) with a noise level a > 0.25, the 
particles first converge to a noisy ring and then switch to 
the rotating state due to the effect of noise. The simulations 
suggest that the transition to the rotating state occurs once 
the velocities of the particles cross an alignment threshold. 
The system, in fact, displays hysteresis: one can force the 
swarm to transition from the ring state with a = 0.2 to the 
rotating state by raising the noise to a = 0.25; however, 
it seems that the inverse transition, i.e. making the swarm 
transition back to the ring state by lowering the noise level, 
is extremely unlikely. 

V. CONCLUSIONS 

To summarize, we studied the dynamics of a self- 
propelling swarm in the presence of noise and a constant 
communication time delay and prove that the delay induces 
a transition that depends upon the size of the interaction 
coupling coefficient. Although our analytical and numerical 
results were obtained using a model with linear, attractive 
interactions, the analysis may be applied to models with more 
general forms of social interaction. 

We further uncovered a complete analytical description of 
the bifurcation point which control the instabilities arising 
from noise induced transitions. The analysis allows us to 
completely classify, using mean field approximations, where 
the swarm exhibits a stable translation, stationary center of 
mass, or rotation. 

In general, our results provide insight into the stability of 
complex systems comprised of individuals interacting with 
one another with a finite time delay in a noisy environment. 
Furthermore, the results may prove to be useful in controlling 
man-made vehicles where actuation and communication are 
delayed, as well as in understanding swarm alignment in 
biological systems. 
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